Quasi-Particle density of states of disordered d-wave superconductors 
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We present a numerical study of the quasi-particle density of states (DoS) of two-dimensional d- 
wave superconductors in the presence of smooth disorder. We find power law scaling of the DoS with 
an exponent depending on the strength of the disorder and the superconducting order parameter 
in quantitative agreement with the theory of Nersesyan et al. (Phys. Rev. Lett. 72, 2628 (1994)). 
For strong disorder a transition to a constant DoS occurs. Our results are in contrast to the case of 
short-ranged disorder. 
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In recent years, motivated by its relevance to the 
physics of the cuprates, d-wave superconductivity has be- 
come a subject of intensive research. The key feature 
distinguishing d-wave superconductors from their s-wave 
relatives is the existence of four zero energy 'nodes' on 
the Fermi surface, in the vicinity of which low energy 
quasi-particle excitations exist. This fermion system — 
four species of relativistic quasi-particles, subject to weak 
static disorder — has vital influence on low energy trans- 
port and thermodynamic properties and, therefore, must 
be a central element of any comprehensive theory of the 
d-wave superconductor. 

Irritatingly, it has proven excruciatingly difficult to 
reach consensus on even basic characteristics of these 
states. The extent of disagreement is most clearly dis- 
played in the debate on the energy dependent mean 
quasi-particle density of states (DoS), p{E): On the 
one hand, application of self-consistent approximation 
schemes or non-perturbative approaches specific to 
certain realizations of the disorder [^-Q, has led to the 
prediction of a finite or even diverging quasi-particle DoS 
at zero energy. In contrast, field theory approaches to 
the problem ||||] categorically predict p(E) = 0, as in 
the non-disordered Dirac system. Yet, even these theo- 
ries within themselves come to varying conclusions as to 
the energy dependence of the DoS for E ^ 0. 

Recently it has become clear that much of this 
controversy roots in the fact that, unlike with more 
conventional disordered fermion systems, the standard 
paradigm of 'insensitivity of global observables to mi- 
croscopic details of the disorder' is apparently violated 
in the d-wave system. Broadly speaking, two different 
categories of disorder have to be distinguished: (i) hard 
scattering off s-wave impurities, mixing the formerly iso- 
lated four low energy quasi-particle sectors, and (ii) soft 
scattering which predominantly leads to randomisation 
within these sectors. Which of these categories is more 
relevant to the physics of the cuprates is not straight- 
forward to decide (see, however, our comments below), 
and both have been investigated theoretically. As for (i), 
there is now overwhelming evidence that, apart from the 



case of asymptotically strong impurities (impurities at 
the 'unitary limit') ||,|,@@, the DoS vanishes as E -> 
|B||Jl3] . Away from zero energy a variety of different DoS 
profiles, depending on the realization of the disorder dis- 
tribution, exist. 

The subject of this Letter is a numerical study of the 
complementary case, (ii). What makes this regime spe- 
cial, and why is it necessary to discriminate between (i) 
and (ii) at all? The distinguishing feature of the soft 
scattering regime is that the four low energy sectors are 
decoupled. The absence of inter-node coupling has pro- 
found and qualitative influence on the low energy prop- 
erties of the system. Indeed |3] , it is the nodal coupling 
criterion, and not so much the specifics of a short range 
correlated disorder distribution jfl| , that holds responsi- 
ble for much of the discrepancy between the field theory 
approaches to the problem. 

Before turning to our numerical analysis of the quasi- 
particle spectrum for soft scattering, let us briefly sum- 
marize some key features of the system: Each node ac- 
commodates a system of Dirac fermions subject to a 
random vector potential which describes the stochastic 
low momentum transfer scattering. For low energies and 
finite size systems - the 'zero-dimensional' limit - the 
properties of the system become fully universal and can 
be described in terms of a suitably constructed random 
matrix theory (RMT) |1| (a). Due to the non-standard 
symmetries of random gauge Dirac fermions (symmetries 
of class Alll in the terminology of Rcf. |lif| ) this theory 
differs profoundly from standard Wigner-Dyson RMT. 
The opposite extreme, thermodynamically extended sys- 
tems is widely accessible to analytical approaches, too 
p5t^6| . Specifically, for the d-wave problem Nersesyan, 
Tsvelik and Wenger (NTW) [§ have shown that the DoS 
scales as (b) 
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where W is the strength of the disorder, and t and A are 
the tight binding coupling strength and order parameter 
of the superconductor, respectively. It has been argued 
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jlTj that at g — 1, i.e. at the zero of the exponent in ([!]), 
a transition to a qualitatively different phase takes place. 
On the strong disorder side of this transition, g > 1, the 
DoS is expected to be energy independent, p(E) = const, 
(c). Further, (d), any amount of hard scattering coupling 
the low energy nodes represents a marginally relevant 
perturbation |l2| driving the system towards the coupled 
regimes (i) mentioned above. 

For completeness we mention that all these features 
find their common origin in the fact that the low energy 
physics of isolated nodes is described by a Wess-Zumino- 
Witten (WZW) model on a group manifold that depends 
on the treatment of the disorder (replica @,Q or super- 
symmetry 0]). In NTW's analysis of this connection 
it has erroneously been assumed that the WZW action 
globally describes the low-energy quasi-particle system, 
independent of the form factor of the scattering. In fact, 
however, the WZW model is readily destabilized by inter- 
node scattering which is one way of explaining the afore- 
mentioned qualitative differences between cases (i) and 
(ii). For a detailed account of the WZW- formulation, 
and its destruction, we refer to Ref. [Q . 

Below we will put the phenomenology (a-d) to a nu- 
merical test. Before turning to a more detailed descrip- 
tion of our analysis, let us summarize the main results. 
We find that the large scale structure of the energy de- 
pendent DoS can be characterized in terms of three dif- 
ferent regimes: For low energies E above the chemical 
potential (owing to the Dirac structure of the problem, 
the DoS is symmetric p(E) = p(—E)) the DoS profile 
is dominated by finite size effects. At E = the DoS 
vanishes in a way which (for strong enough disorder) is 
described by RMT. The extent of the low energy regime 
shrinks with increasing system size. For larger energies, 
it is succeeded by a regime of power law scaling. Vary- 
ing the two basic parameters A/t and W/t characterizing 
the model, we find agreement with eq. (Q). For disorder 
strength in excess of g = 1, the DoS assumes a constant 
value, in accord with the prediction of Ref. O]. The 
scaling regime ends at energies of the order of A where a 
non-universal high energy regime, not considered in this 
Letter, begins. Upon lowering the correlation length of 
the disorder, the scaling behaviour observed in the center 
portions of the band is rapidly destructed. 

We consider the lattice quasi-particle Hamiltonian 
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with the hopping matrix elements iy , chemical potential 
/i, and order parameter Ajj. The sums run over points of 
a two-dimensional square lattice with spacing a and the 
operators c\ a create a spin-1/2 particle of spin a at site i. 
In the following, we take only into account on-site poten- 
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FIG. 1. Density of states for disorder of strength 
W — 0,1,..., 8, 10 (bottom to top at E = 0), correlation 
length £ = 2, and A = 1. The system size is L = 33 and 
the broadening T = 0.05. The inset shows the same data on 
a smaller scale with F = 0.0005. The finite DoS at E = is 
due to the finite broadening I\ 

nience we set p = (the half-filled band.) The order pa- 
rameter Ay = A(5i t j± eal —Si_j± £y ) has d x 2_ y 2 -symmetry. 
In the following we will only consider disorder in the on- 
site potentials and not determine the order parameter 
self-consistently. (As pointed out in jllj , a self-consistent 
determination of the order parameter may be necessary 
to quantitatively compare with experimental data.) Fi- 
nite correlations in the disorder potential are introduced 
through 
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where is the unit vector in /c-direction. For convc- 



£ = cx P(~2|rj | 2 /£ 2 ), where the fj are independent 
random variables, uniformly distributed in the interval 
[—1/2, 1/2], and IT is a measure of the strength of the 
disorder. Taking £ = 2a the coupling of the four Dirac 
nodes is much weaker than the coupling of states in the 
vicinity of a single node (by a factor of exp(— 2ir 2 )). 

In the following, we take the lattice constant a and 
the hopping amplitude t as the measures of length and 
energy, respectively. The linear dimensions of the system 
vary between 15 and 45 and the spectrum is averaged 
over 64 points in the first Brioullin zone. We diagonalize 
the Hamiltonian (& for 16 to 56 disorder realizations 
and calculate the DoS p(E) = L~ 2 JT Sr(E — Ei), where 
Sr(E) is a normalized Gaussian of width T. A finite T 
smoothes the DoS, but also washes out narrow features. 

Fig. [I] shows the DoS for various values of disorder. 
The three regimes mentioned earlier are the low-energy 
region where the bump develops, the intermediate regime 
up to the approximate crossing point near E = 1, and 
the high-energy regime beyond. We will first discuss the 
most interesting intermediate regime before turning our 
attention to low energies. In order to compare our data to 
the results of NTW, we fit the data to power laws p(E) oc 
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FIG. 2. Double logarithmic plot of the density of states of 
Fig. |l|. Disorder ranges from W = 1 to 10. Dots (•) repre- 
sent data and lines power law fits to the respective intervals. 
Inset: Density of states for W = 2 and L = 15 (dotted), 25 
(short-dashed), 35 (long-dashed), and 45 (solid). Note that 
the numerical uncertainties are considerably smaller than the 
amplitude of the fluctuations. 
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FIG. 3. Exponents a extracted from the fitted curves in 
Fig. |as a function of disorder W for A = 1. The solid curve 
is the result of NTW, eq. <JTJ> . 

E a in a interval [E m i n ,E max \. £Wn is cfioosen such as 
to exclude the first maximum. From the inset in Fig. ^ 
it is clear that this feature as well as the fluctuations 
are finite-size effects that appear to vanish in the large 
system limit. The upper limit E max is a high energy cut- 
off of the order of A beyond which the NTW theory no 
longer applies. At the systems sizes considered in this 
work, the ratio E max / E min is about 5. This is a rather 
narrow range to establish a power law. Nevertheless, we 
feel that our procedure is justified in the present case, as 
it is not just a single power law with a single exponent 
that we are dealing with. Instead, eq. ([!]) predicts that 
there is a whole family of power laws with the exponents 
depending in a unique way on the parameters W and A. 
It is the agreement of this whole functional dependence 
that gives us confidence in the validity of our analysis 
even when the establishment of every single power law 
might be questionable. Figure || shows the exponents a 
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FIG. 4. Density of states for order parameter A = 0.1 to 
1.0 (top to bottom) and disorder W = 3. Each curve is shifted 
by a factor of 1.2 for clarity. 

extracted from the fitted curves in Fig. || together with 
the result of eq. ([!]) of || . A good agreement is apparent 
up to a disorder strength of W w 7 where the NTW 
exponent changes sign. Numerically, we do not find a 
divergent DoS at stronger disorder but rather a finite 
value (a — 0) as predicted by Gurarie [p7| . 

To further test eq. (0) we fix the disorder strength W — 
3 and vary the order parameter A. Figure ^ shows the 
DoS for A between 0.1 and 1.0 as a function of E/A. 
This rescaling takes care of the fact that the mean level 
spacing of the clean system is proportional to A. The 
anisotropy dependence of the fitted exponents is shown 
in Fig. ||. Again, reasonably good agreement with eq. ([!]) 
is found for A > 0.2 while a constant and not a diverging 
DoS is found at A = 0.1. 

The inset of Fig. [l] shows that the DoS does indeed 
vanish at E = 0. At weak disorder (W < 4) we see a 
remnant of the clean spectrum. Here the disorder is too 
weak to couple neighboring states in momentum space. 
At W ss 5 this coupling exceeds the level separation and 
the universal RMT behaviour for systems of class AIII, 
with a DoS 'microgap' linear in energy, develops ||l9|| . 

To conclude the central part of our analysis, we have 
presented a numerical study of the quasi-particle spec- 
trum in d-wave superconductors with soft scattering. In 
the regime of moderate disorder strength, we obtain scal- 
ing of the DoS that agrees quantitatively with the ana- 
lytical results of NTW. For strong disorder, the analy- 
sis confirms the recent prediction [jlTj of an energetically 
constant background DoS. In a way our analysis is com- 
plementary to recent studies of hard scattering systems 
Jll| and we are left with the question which of these alter- 
natives is more relevant to the physics of 'real' cuprates. 
It is probably difficult to give a universally applicable 
answer to this question. On the face of it, disorder in 
high T c superconductors is due to small metallic donors, 
e.g. Zn- impurities, in favour of the hard variant, (i). 
(Concrete evidence for the presence of hard scattering in 
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FIG. 5. Exponents a for W = 3 as a function of order 
parameter A. The solid curve is the result of NTW, eq. ([[]). 

Zn-doped systems is provided by the experimental ob- 
servation |2(J of so-called mid-gap resonances [i.e. res- 
onances due to bound states forming in the immediate 
vicinity of a strong local impurity]). On the other hand, 
the very existence of a <i-wave phase, stabilized through 
a mechanism unknown at present, would not be compat- 
ible with too strong an amount of — pairbreaking — hard 
scattering, i.e. the renormalized effective potential seen 
by the quasi-particle states may well carry characteristics 
of type (ii) and be predominantly soft (see Refs. J2l[ for 
a more elaborate discussion of this point.) Equally im- 
portant, the net features of the quasi-particle system are 
not only determined by the fixed microscopic structure of 
the disorder background but also depend on temperature 
and observation energy. E.g. consider a system with a 
certain residual amount of hard scattering superimposed 
on a predominantly soft background. For quasi-particle 
energies in excess of the the inter-node scattering rate, 
the coupling between the nodes is inessential and the 
characteristics of the soft system will prevail. Lowering 
the energies, a crossover towards the hard system takes 
place. 

Although the picture above suggests, that 'real life' 
systems will typically display complex crossover be- 
haviour, recent progress in purely theoretical under- 
standing of disorder in d-wave superconductors has been 
tremendous. It seems to be clear now, that much of 
the controversy that developed around the profile of the 
quasi-particle DoS is related to non-congruent modellings 
of the disorder. Indeed, the majority of theoretical ap- 
proaches to the problem sits comfortably with one of 
the disorder realizations investigated numerically in this 
work or in complementary papers [ j^jOj . (There is one 
prominent exception to that rule, viz. Rcfs. || where 
a finite DoS for continuously distributed disorder was 
predicted. To our understanding the discrepancy is ex- 
plained by the peculiar lattice implementation underlying 
these papers (see Ref. Jl2| for a more elaborate discus- 
sion of this point) which implies that no superconductor 



is modelled.) Broadly speaking, there seem to be three 
categories of disorder that have to be distinguished: bi- 
nary alloy type scatterer at the unitary limit, large mo- 
mentum transfer scatterers of non-unitary type and soft 
disorder. In spite of the relative — as compared to nor- 
mal metals — complexity of this classification, the cate- 
gories within themselves still display a large amount of 
universality. 

Acknowledgement: We than B.D. Simons for discus- 
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